trimmedAMKE_probability of nest failure

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(knitr)
library(jagsUI)

## Warning: package 'jagsUI' was built under R version 4.4.3

library(DiagrammeR)

## Warning: package 'DiagrammeR' was built under R version 4.4.3

set.seed(1019)

Seed is set to 1019 for reproducibility.

Data

Manual data, pen-wise, inclusive of multiple broods.

dflongforcount <- data.frame(
  dose=factor(c(rep("cont",6),rep("low",7),rep("med",6),rep("high",8)),
              levels=c("cont","low","med","high")),
  dosenum=c(rep(0.000121,6),rep(0.121,7),rep(0.45,6),rep(2.04,8)),
  pen=c(31,32,33,34,35,35,36,37,38,38,39,40,40,41,42,43,43,44,45,46,46,46,47,48,48,48,50),
  clutch=c(1,1,1,1,1,2,1,1,1,2,1,1,2,1,1,1,2,1,1,1,2,3,1,1,2,3,1),
  numhatched=c(3,2,2,1,0,0,4,1,0,3,3,0,0,3,0,0,0,1,1,0,0,0,2,0,0,0,2)
)

dflongforcount$failure <- ifelse(dflongforcount$numhatched==0,1,0)

JAGS setup and run of Poisson GLM

jags.data <- list(y=dflongforcount$numhatched,
                  conc=dflongforcount$dosenum,
                  n=nrow(dflongforcount))

cat(file="model_AMKE_poisson_eggs hatched by normal dose.txt","
  model{
  #priors and linear models
  alpha ~ dnorm(log(3),100) # mean and precision in jags precision is inverse of the variance, so 1/sigma^2 and in this case suggests the priors are mean=0, sd=0.1 (exp(0.1) on the real scale) (sqrt(1/100))
  beta ~ dnorm(0,0.001)
    
  #likelihood for the poisson GLM
  for(i in 1:n){
    y[i] ~ dpois(lambda[i]) #stochastic
    log(lambda[i]) <- alpha + beta * conc[i]
  }
  mean.exp.count <- exp(alpha)
}
")

inits <- function(){list(alpha=rnorm(1,0,1),beta=rnorm(1,0,1))}

parameters <- c("alpha","beta","mean.exp.count","lambda")

ni <- 50000
nb <- 10000
nc <- 3
nt <- 10
na <- 1000

out1 <- jags(jags.data, inits, parameters, "model_AMKE_poisson_eggs hatched by normal dose.txt", n.iter=ni, n.burnin=nb, n.chains=nc, n.thin=nt, n.adapt=na, parallel=T)

## 
## Processing function input....... 
## 
## Done. 
##  
## Beginning parallel processing using 3 cores. Console output will be suppressed.
## 
## Parallel processing completed.
## 
## Calculating statistics....... 
## 
## Done.

test <- data.frame(out1$sims.list[1:2])
#write_csv(test, "out1posterioralphaandbeta.csv")

\[ \alpha \sim N(3,0.1)\\ \beta \sim N(0,31.62)\\ y \sim Poisson(\lambda)\\ log(\lambda) = alpha + \beta \sum PFAS \]

JAGS output wrangling for Poisson GLM

contm <- out1$mean$lambda[1] #2.575
lowm <- out1$mean$lambda[7] #2.264
medm <- out1$mean$lambda[14] #1.605
highm <- out1$mean$lambda[20] #0.343

contl <- out1$q2.5$lambda[1] #2.164
lowl <- out1$q2.5$lambda[7] #1.906
medl <- out1$q2.5$lambda[14] #1.181
highl <- out1$q2.5$lambda[20] #0.081

contu <- out1$q97.5$lambda[1] #3.041
lowu <- out1$q97.5$lambda[7] #2.66
medu <- out1$q97.5$lambda[14] #2.034
highu <- out1$q97.5$lambda[20] #0.812


### probability of having a certain number of viable eggs

counteggdf <- data.frame(
  dose=c(0.000121, 0.121, 0.45, 2.04),
  mean=c(contm, lowm, medm, highm),
  lower=c(contl, lowl, medl, highl),
  upper=c(contu, lowu, medu, highu)
)

effectzerodf <- data.frame(
  dose=c(0.000121, 0.121, 0.45, 2.04),
  mean=c(dpois(0, lambda=contm),dpois(0, lambda=lowm),dpois(0, lambda=medm),dpois(0, lambda=highm)),
  lower=c(dpois(0, lambda=contl),dpois(0, lambda=lowl),dpois(0, lambda=medl),dpois(0, lambda=highl)),
  upper=c(dpois(0, lambda=contu),dpois(0, lambda=lowu),dpois(0, lambda=medu),dpois(0, lambda=highu)),
  hatchedeggs=rep(0,4)
)

effect3df <- data.frame(
  dose=c(0.000121, 0.121, 0.45, 2.04),
  mean=c(dpois(3, lambda=contm),dpois(3, lambda=lowm),dpois(3, lambda=medm),dpois(3, lambda=highm)),
  lower=c(dpois(3, lambda=contl),dpois(3, lambda=lowl),dpois(3, lambda=medl),dpois(3, lambda=highl)),
  upper=c(dpois(3, lambda=contu),dpois(3, lambda=lowu),dpois(3, lambda=medu),dpois(3, lambda=highu)),
  hatchedeggs=rep(3,4)
)

effect2df <- data.frame(
  dose=c(0.000121, 0.121, 0.45, 2.04),
  mean=c(dpois(2, lambda=contm),dpois(2, lambda=lowm),dpois(2, lambda=medm),dpois(2, lambda=highm)),
  lower=c(dpois(2, lambda=contl),dpois(2, lambda=lowl),dpois(2, lambda=medl),dpois(2, lambda=highl)),
  upper=c(dpois(2, lambda=contu),dpois(2, lambda=lowu),dpois(2, lambda=medu),dpois(2, lambda=highu)),
  hatchedeggs=rep(2,4)
)

effect5df <- data.frame(
  dose=c(0.000121, 0.121, 0.45, 2.04),
  mean=c(dpois(5, lambda=contm),dpois(5, lambda=lowm),dpois(5, lambda=medm),dpois(5, lambda=highm)),
  lower=c(dpois(5, lambda=contl),dpois(5, lambda=lowl),dpois(5, lambda=medl),dpois(5, lambda=highl)),
  upper=c(dpois(5, lambda=contu),dpois(5, lambda=lowu),dpois(5, lambda=medu),dpois(5, lambda=highu)),
  hatchedeggs=rep(5,4)
)

effect1df <- data.frame(
  dose=c(0.000121, 0.121, 0.45, 2.04),
  mean=c(dpois(1, lambda=contm),dpois(1, lambda=lowm),dpois(1, lambda=medm),dpois(1, lambda=highm)),
  lower=c(dpois(1, lambda=contl),dpois(1, lambda=lowl),dpois(1, lambda=medl),dpois(1, lambda=highl)),
  upper=c(dpois(1, lambda=contu),dpois(1, lambda=lowu),dpois(1, lambda=medu),dpois(1, lambda=highu)),
  hatchedeggs=rep(1,4)
)

effect4df <- data.frame(
  dose=c(0.000121, 0.121, 0.45, 2.04),
  mean=c(dpois(4, lambda=contm),dpois(4, lambda=lowm),dpois(4, lambda=medm),dpois(4, lambda=highm)),
  lower=c(dpois(4, lambda=contl),dpois(4, lambda=lowl),dpois(4, lambda=medl),dpois(4, lambda=highl)),
  upper=c(dpois(4, lambda=contu),dpois(4, lambda=lowu),dpois(4, lambda=medu),dpois(4, lambda=highu)),
  hatchedeggs=rep(4,4)
)

bigeggprobdf <- rbind(
  effectzerodf,
  effect1df,
  effect2df,
  effect3df,
  effect4df,
  effect5df
)

Table of Poisson GLM output

Probability of having a certain number of hatched eggs at each dose with confidence intervals.
Note that upper and lower are strange because some dose response are decreasing and some are increasing depending on the number of hatched eggs.

bigeggprobdf |>
  #filter(hatchedeggs==0)|>
  kable(, digits=5)

dose	mean	lower	upper	hatchedeggs
0.00012	0.07615	0.11519	0.04752	0
0.12100	0.10421	0.14819	0.06999	0
0.45000	0.20227	0.30665	0.13041	0
2.04000	0.71318	0.92599	0.44442	0
0.00012	0.19609	0.24894	0.14477	1
0.12100	0.23566	0.28293	0.18614	1
0.45000	0.32326	0.36247	0.26566	1
2.04000	0.24107	0.07120	0.36042	1
0.00012	0.25247	0.26901	0.22053	2
0.12100	0.26645	0.27010	0.24750	2
0.45000	0.25831	0.21423	0.27058	2
2.04000	0.04074	0.00274	0.14615	2
0.00012	0.21671	0.19379	0.22396	3
0.12100	0.20084	0.17190	0.21940	3
0.45000	0.13761	0.08441	0.18373	3
2.04000	0.00459	0.00007	0.03951	3
0.00012	0.13951	0.10470	0.17058	4
0.12100	0.11354	0.08205	0.14586	4
0.45000	0.05498	0.02495	0.09357	4
2.04000	0.00039	0.00000	0.00801	4
0.00012	0.07185	0.04526	0.10394	5
0.12100	0.05135	0.03133	0.07758	5
0.45000	0.01757	0.00590	0.03812	5
2.04000	0.00003	0.00000	0.00130	5

Interpolate adverse increase in probability of nest failure using Poisson GLM

Using the posterior draws, find the doses bracketing 95% probability of observing that amount of nest failure.

Background probability of having zero hatched eggs

dpois(0,lambda=contm)

## [1] 0.07614831

0.2 added risk of having zero hatched eggs

dpois(0,lambda=contm)+0.2

## [1] 0.2761483

Highest concentration where probability of having zero eggs is <0.28

i.e. “safe” concentration

### do prediction interpolation for having zero

alt.conc.pred <- seq(from=0.0001,to=3, length.out=1000)
pred2 <- array(NA,dim=c(length(alt.conc.pred),12000))
for(i in 1:12000){
  pred2[,i] <- exp(out1$sims.list$alpha[i]+out1$sims.list$beta[i]*alt.conc.pred)
}

pm2 <- apply(pred2,1, mean)
pm2l <- apply(pred2,1, quantile, probs=0.975)
pm2u <- apply(pred2,1, quantile, probs=0.025)

prob <- apply(pred2, 1, function(x) mean((x-1)<0))

dprob0 <- apply(pred2,1,function(x) mean(dpois(0,lambda=x)))
dprob0min <- apply(pred2,1,function(x) quantile(dpois(0,lambda=x), probs=0.05))
dprob0max <- apply(pred2,1,function(x) quantile(dpois(0,lambda=x), probs=0.95))

#plot(dprob0~alt.conc.pred, ylim=c(0,1))
#points(dprob0min~alt.conc.pred, col="darkgray")
#points(dprob0max~alt.conc.pred, col="darkgray")


mintwentythreshold <- max(alt.conc.pred[dprob0min<0.28])
maxtwentythreshold <- max(alt.conc.pred[dprob0max<0.28])
meantwentythreshold <- max(alt.conc.pred[dprob0<0.28])
mintwentythreshold

## [1] 1.090154

meantwentythreshold

## [1] 0.648727

maxtwentythreshold

## [1] 0.432518

Plot of predicted number of eggs

par(mai=c(2,1.2,0.1,0.3))
plot(jitter(numhatched, 0.5)~jitter(as.numeric(factor(dose))), data=dflongforcount, pch=3, ylab="Number of hatched eggs", xlab=expression(Adult~daily~dose~mg/kg~Sigma*PFAS),
     bty="n", las=1, axes=F, cex.lab=1.5)
axis(1, at=c(1,2,3,4), labels=c("0.000121", "0.121", "0.45", "2.04"), font=1, cex.axis=1.5)
axis(2, at=c(0,1,2,3,4), labels=c(0,1,2,3,4), font=1, las=1, cex.axis=1.5)
arrows(as.numeric(factor(counteggdf$dose)),counteggdf$lower,as.numeric(factor(counteggdf$dose)),counteggdf$upper, lwd=2, length=0.1, angle=90, code=3, col=c("darkgreen","blue","darkorange","red"))
points(mean~as.numeric(factor(dose)), data=counteggdf, type="l", col="darkgray", lwd=2)
points(mean~as.numeric(factor(dose)), data=counteggdf, type="p", col=c("darkgreen","blue","darkorange","red"), pch=16, cex=1.5)
points(c(1.1,1.1),c(2.95,3.05),col="purple", pch=2)
axis(1, at=c(1,2,3,4), labels=c("0.003", "3.35", "8.34", "65.84"), font=1, cex.axis=1.5, line=5)
mtext(expression(Mean~egg~mg/kg~Sigma*PFAS), side=1, line=8, font=1, cex=1.5)

par(mai=c(1,1,0.1,0.1))
plot(pm2~alt.conc.pred, log="", ylim=c(0,3), type="l", bty="n", xlim=c(-0.25,3),
     ylab="Mean number of hatched eggs", xlab=expression(paste("Adult daily dose ", Sigma,"PFAS mg/kg")), las=1,
     font=1, font.lab=1, pch=16, cex=2, lwd=2, cex.lab=1.5, cex.axis=1.5)
points(pm2l~alt.conc.pred, type="l", lty=2, col="darkgray", lwd=2)
points(pm2u~alt.conc.pred, type="l", lty=2, col="darkgray", lwd=2)
text(-0.35,2.575,"2.575", pos=4, col="darkgreen")
text(-0.35,1.264456,"1.264", pos=4, col="red")
segments(0,1.264456,1.5,1.264456, col="red")
segments(0.6757532,1.264456,0.6757532,0.5, col="black")
text(0.6757532,0.5,"0.68 mg/kg", pos=1, col="black", xpd=NA)
segments(0.4054919,1.264456,0.4054919,1.064456, col="black")
#text(0.432518,0.48,"0.43 mg/kg\nP(0 eggs)<0.28", pos=3, col="black")
text(0.4054919,1.064456,"0.41 mg/kg", pos=1, col="black")
segments(1.237296,1.264456,1.237296,1.464456, col="black")
#text(1.090154,0.18,"1.09 mg/kg\nP(0 eggs)<0.28", pos=1, col="black", xpd=NA)
text(1.237296,1.464456,"1.24 mg/kg", pos=3, col="black", xpd=NA)

max(alt.conc.pred[pm2>1.264456])

## [1] 0.6757532

max(alt.conc.pred[pm2l>1.264456])

## [1] 1.237296

max(alt.conc.pred[pm2u>1.264456])

## [1] 0.4054919

par(mai=c(1.2,1.2,0.1,0.75))
plot(dprob0~alt.conc.pred, log="", ylim=c(0,1), type="l", bty="n", xlim=c(-0.25,3),
     ylab="Probability of having 0 hatched eggs", xlab=expression(paste("Adult daily dose ", Sigma,"PFAS mg/kg")), las=1,
     font=1, font.lab=1, pch=16, cex=2, lwd=2, cex.lab=1.5, cex.axis=1.5)
points(dprob0min~alt.conc.pred, type="l", lty=2, col="darkgray", lwd=2)
points(dprob0max~alt.conc.pred, type="l", lty=2, col="darkgray", lwd=2)
abline(h=c(0,1))
text(-0.35,dpois(0,lambda=contm),"0.08", pos=4, col="darkgreen")
text(-0.35,dpois(0,lambda=contm)+0.2,"0.28", pos=4, col="red")
segments(0,dpois(0,lambda=contm)+0.2,1.5,dpois(0,lambda=contm)+0.2, col="red")
segments(0.432518,0.28,0.432518,0.48, col="black")
#text(0.432518,0.48,"0.43 mg/kg\nP(0 eggs)<0.28", pos=3, col="black")
text(0.432518,0.48,"0.43 mg/kg", pos=3, col="black")
segments(1.090154,0.28,1.090154,0.18, col="black")
#text(1.090154,0.18,"1.09 mg/kg\nP(0 eggs)<0.28", pos=1, col="black", xpd=NA)
text(1.090154,0.18,"1.09 mg/kg", pos=1, col="black", xpd=NA)
segments(0.648727,0.28,0.648727,0.08, col="black")
text(0.648727,0.08,"0.65 mg/kg", pos=1, col="black", xpd=NA)

JAGS setup and run of Bernoulli GLM

jags.data2 <- list(y=dflongforcount$failure,
                  conc=dflongforcount$dosenum,
                  n=nrow(dflongforcount))

cat(file="model_AMKE_bernoulli_eggs hatched by normal dose.txt","
  model{
  #priors and linear models
  alpha <- logit(mean.theta)
  mean.theta ~ dunif(0,0.5)
  beta ~ dnorm(0,0.001)
    
  #likelihood for the bernoulli GLM
  for(i in 1: n){
    y[i] ~ dbern(theta[i]) #stochastic
    logit(theta[i]) <- alpha + beta * conc[i]
  }
}
")

inits2 <- function(){list(mean.theta=runif(1,0, 0.5),beta=rnorm(1,0,1))}

parameters2 <- c("alpha","beta","mean.theta","theta")

ni <- 50000
nb <- 10000
nc <- 3
nt <- 10
na <- 1000

out2 <- jags(jags.data2, inits2, parameters2, "model_AMKE_bernoulli_eggs hatched by normal dose.txt", n.iter=ni, n.burnin=nb, n.chains=nc, n.thin=nt, n.adapt=na, parallel=T)

## 
## Processing function input....... 
## 
## Done. 
##  
## Beginning parallel processing using 3 cores. Console output will be suppressed.
## 
## Parallel processing completed.
## 
## Calculating statistics....... 
## 
## Done.

#jagsUI::traceplot(out2)

test2 <- data.frame(out2$sims.list[1:2])
#write_csv(test, "out1posterioralphaandbeta.csv")

\[ \alpha \sim \text{Uniform}(0,0.5)\\ \beta \sim N(0,31.62)\\ y \sim \text{Bernoulli}(\theta)\\ logit(\theta) = alpha + \beta \sum PFAS \]

JAGS output wrangling for Bernoulli GLM

contb <- out2$mean$theta[1]
lowb <- out2$mean$theta[7]
medb <- out2$mean$theta[14]
highb <- out2$mean$theta[21]

contbl <- out2$q2.5$theta[1]
lowbl <- out2$q2.5$theta[7]
medbl <- out2$q2.5$theta[14]
highbl <- out2$q2.5$theta[21]

contbu <- out2$q97.5$theta[1]
lowbu <- out2$q97.5$theta[7]
medbu <- out2$q97.5$theta[14]
highbu <- out2$q97.5$theta[21]

### probability of a failed nest (having zero hatched eggs)

failuredf <- data.frame(
  dose=c(0.000121, 0.121, 0.45, 2.04),
  mean=c(contb, lowb, medb, highb),
  lower=c(contbl, lowbl, medbl, highbl),
  upper=c(contbu, lowbu, medbu, highbu)
)



# plot(jitter(dflongforcount$failure, 0.2)~jitter(as.numeric(factor(dflongforcount$dosenum)),0.5), axes=F, pch=3, ylab="Per nest probability of failure", xlab="Diet concentration (mg/kg-d sumPFAS)", )
# lines(as.numeric(factor(dflongforcount$dosenum)), out2$mean$theta, col="blue", lwd=3)
# lines(as.numeric(factor(dflongforcount$dosenum)), out2$q2.5$theta, lwd=3, lty=2)
# lines(as.numeric(factor(dflongforcount$dosenum)), out2$q97.5$theta, lwd=3, lty=2)
# axis(1,at=unique(as.numeric(factor(dflongforcount$dosenum))), labels=sort(unique(round(dflongforcount$dosenum,2))))
# axis(2,las=1)

Plot of predicted probability of failure (zero hatched eggs)

par(mai=c(2,1.2,0.1,0.3))
plot(jitter(failure, 0.5)~jitter(as.numeric(factor(dose))), data=dflongforcount, pch=3, ylab="Probability of nest failure", xlab=expression(Adult~daily~dose~mg/kg~Sigma*PFAS),
     bty="n", las=1, axes=F, cex.lab=1.5)
axis(1, at=c(1,2,3,4), labels=c("0.000121", "0.121", "0.45", "2.04"), font=1, cex.axis=1.5)
axis(2, at=c(0,0.5,1), labels=c(0,0.5,1), font=1, las=1, cex.axis=1.5)
arrows(as.numeric(factor(failuredf$dose)),failuredf$lower,as.numeric(factor(failuredf$dose)),failuredf$upper, lwd=2, length=0.1, angle=90, code=3, col=c("darkgreen","blue","darkorange","red"))
points(mean~as.numeric(factor(dose)), data=failuredf, type="l", col="darkgray", lwd=2)
points(mean~as.numeric(factor(dose)), data=failuredf, type="p", col=c("darkgreen","blue","darkorange","red"), pch=16, cex=1.5)
points(c(0.95,1.05),c(0,0),col="purple", pch=2)
axis(1, at=c(1,2,3,4), labels=c("0.003", "3.35", "8.34", "65.84"), font=1, cex.axis=1.5, line=5)
mtext(expression(Mean~egg~mg/kg~Sigma*PFAS), side=1, line=8, font=1, cex=1.5)

Interpolate adverse increase in probability of nest failure using Bernoulli GLM

Using the posterior draws, find the doses bracketing 95% probability of observing that amount of nest failure.

Background probability of having zero hatched eggs

dbinom(1,1,prob=contb)

## [1] 0.3577885

0.2 added risk of having zero hatched eggs

dbinom(1,1,prob=contb)+0.2

## [1] 0.5577885

Highest concentration where probability of having zero eggs is <0.28

i.e. “safe” concentration

### do prediction interpolation for having zero

alt.conc.pred <- seq(from=0.0001,to=3, length.out=1000)
pred3 <- array(NA,dim=c(length(alt.conc.pred),12000))
for(i in 1:12000){
  pred3[,i] <- plogis(out2$sims.list$alpha[i]+out2$sims.list$beta[i]*alt.conc.pred)
}

pm3 <- apply(pred3,1, mean)
pm3l <- apply(pred3,1, quantile, probs=0.975)
pm3u <- apply(pred3,1, quantile, probs=0.025)

prob2 <- apply(pred3, 1, function(x) mean((x-1)<0))

dprob02 <- apply(pred3,1,function(x) mean(dbinom(1,1,prob=x)))
dprob02min <- apply(pred3,1,function(x) quantile(dbinom(1,1,prob=x), probs=0.05))
dprob02max <- apply(pred3,1,function(x) quantile(dbinom(1,1,prob=x), probs=0.95))

#plot(dprob0~alt.conc.pred, ylim=c(0,1))
#points(dprob0min~alt.conc.pred, col="darkgray")
#points(dprob0max~alt.conc.pred, col="darkgray")


mintwentythreshold2 <- max(alt.conc.pred[dprob02min<0.5579766])
maxtwentythreshold2 <- max(alt.conc.pred[dprob02max<0.5579766])
meantwentythreshold2 <- max(alt.conc.pred[dprob02<0.5579766])
mintwentythreshold2

## [1] 2.915919

meantwentythreshold2

## [1] 0.8919622

maxtwentythreshold2

## [1] 0.3574454

par(mai=c(1.2,1.2,0.1,0.75))
plot(dprob02~alt.conc.pred, log="", ylim=c(0,1), type="l", bty="n", xlim=c(-0.25,3),
     ylab="Probability of nest failure", xlab=expression(paste("Adult daily dose ", Sigma,"PFAS mg/kg")), las=1,
     font=1, font.lab=1, pch=16, cex=2, lwd=2, cex.lab=1.5, cex.axis=1.5)
points(dprob02min~alt.conc.pred, type="l", lty=2, col="darkgray", lwd=2)
points(dprob02max~alt.conc.pred, type="l", lty=2, col="darkgray", lwd=2)
abline(h=c(0,1))
text(-0.35,dbinom(1,1,prob=contb),"0.36", pos=4, col="darkgreen")
text(-0.35,dbinom(1,1,prob=contb)+0.2,"0.56", pos=4, col="red")
segments(0,dbinom(1,1,prob=contb)+0.2,3,dbinom(1,1,prob=contb)+0.2, col="red")
segments(0.3514396,0.56,0.3514396,0.66, col="black")
text(0.3514396,0.66,"0.35 mg/kg", pos=3, col="black")
segments(3,0.56,3,0.26, col="black")
text(2.915919,0.26,"2.92 mg/kg", pos=1, col="black", xpd=NA)
segments(0.8919622,0.56,0.8919622,0.36, col="black")
text(0.8919622,0.36,"0.89 mg/kg", pos=1, col="black", xpd=NA)

Stochastic population trajectories

Simple model diagram

mermaid("
  graph RL
    A(Adult);
    J(Juvenile);
    I(Immigration);
    
    A--Fecundity-->J;
    A--Adult Survival-->A;
    J--Juvenile Survival-->A;
    I-->A;
", 
width=500, height=500)

Build distributions for demographic rates

Survival and immigration rates from McClure et al. 2021. McClure CJW et al. 2021. Demography of a widespread raptor across disparate regions. Ibis. 163(2):658–670 https://doi.org/10.1111/ibi.12916

Juvenile survival

juvsurvdist <- c(
  rnorm(1000, mean=0.13, sd=0.06),
  rnorm(1000, mean=0.15, sd=0.05),
  rnorm(1000, mean=0.09, sd=0.03),
  rnorm(1000, mean=0.08, sd=0.06)
)
juvsurvdist <- juvsurvdist[!juvsurvdist<0]
par(mfrow=c(1,2),mai=c(0.8,0.8,0.1,0.1))
hist(juvsurvdist, main="")
plot(density(juvsurvdist), main="",bty="n")

par(mfrow=c(1,1),mai=c(1,1,0.5,0.5))
summary(juvsurvdist)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0003585 0.0756261 0.1106364 0.1159136 0.1525845 0.3493614

sd(juvsurvdist)

## [1] 0.05639346

Juvenile survival variability

Across the populations.

sd(c(0.13, 0.15, 0.09, 0.08))

## [1] 0.03304038

Adult survival

adtsurvdist <- c(
  rnorm(1000, mean=0.50, sd=0.03),
  rnorm(1000, mean=0.48, sd=0.03),
  rnorm(1000, mean=0.55, sd=0.03),
  rnorm(1000, mean=0.56, sd=0.08)
)
#adtsurvdist <- adtsurvdist[adtsurvdist>0]
par(mfrow=c(1,2),mai=c(0.8,0.8,0.1,0.1))
hist(adtsurvdist, main="")
plot(density(adtsurvdist), main="",bty="n")

par(mfrow=c(1,1),mai=c(1,1,0.5,0.5))
summary(adtsurvdist)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.3272  0.4826  0.5151  0.5229  0.5548  0.7859

sd(adtsurvdist)

## [1] 0.05953504

Adult survival variability

sd(c(0.50, 0.48, 0.55, 0.56))

## [1] 0.0386221

Immigration

immidist <- c(
  rnorm(1000, mean=0.21, sd=0.06),
  rnorm(1000, mean=0.21, sd=0.05),
  rnorm(1000, mean=0.15, sd=0.04),
  rnorm(1000, mean=0.20, sd=0.09)
)
immidist <- immidist[!immidist<0]
par(mfrow=c(1,2),mai=c(0.8,0.8,0.1,0.1))
hist(immidist, main="")
plot(density(immidist), main="",bty="n")

par(mfrow=c(1,1),mai=c(1,1,0.5,0.5))
summary(immidist)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.001566 0.148561 0.191600 0.194103 0.236840 0.474021

sd(immidist)

## [1] 0.06687488

Immigration variability

sd(c(0.21, 0.21, 0.15, 0.20))

## [1] 0.02872281

Fecundity

fecunddist <- c(
  rnorm(1000, mean=1.55, sd=0.13),
  rnorm(1000, mean=1.36, sd=0.05),
  rnorm(1000, mean=1.38, sd=0.06),
  rnorm(1000, mean=1.32, sd=0.13)
)
#fecunddist <- fecunddist[!fecunddist<0]
par(mfrow=c(1,2),mai=c(0.8,0.8,0.1,0.1))
hist(fecunddist, main="")
plot(density(fecunddist), main="",bty="n")

par(mfrow=c(1,1),mai=c(1,1,0.5,0.5))
summary(fecunddist)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.9304  1.3272  1.3820  1.4032  1.4625  1.9546

sd(fecunddist)

## [1] 0.1328875

Fecundity variability

sd(c(1.55, 1.36, 1.38, 1.32))

## [1] 0.1014479

trial run using Chapter.Section 3.3 from IPM book

3.3 Classical analysis of a matrix population model

3.3.

3.3.5 Analysis of a matrix population model with different sources of stochasticity and parameter uncertainty

See https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R

Code has been modified so “temporal” stochasticity is replaced with “population” stochasticity based on the 4 populations in Mcclure et al. 2021.

R stuff to make it happen, borrowed from IPM book package

# https://github.com/kenkellner/IPMbook/blob/main/R/stopifnot.R

stopifnotProbability <- function(arg, allowNA=FALSE) {
  name <- deparse(substitute(arg))
  if(allowNA && all(is.na(arg))) {  # An all-NA vector is logical, but ok.
    # do nothing
  } else {
    if(!allowNA && any(is.na(arg)))
      stop("Argument '", name, "' must not contain NA or NaN.", call.=FALSE)
    if(!is.numeric(arg))
      stop("Argument '", name, "' must be numeric.", call.=FALSE)
    if(any(arg < 0 | arg > 1, na.rm=TRUE)) {
      if(allowNA) {
        stop("Argument '", name, "' must be a probability between 0 and 1, or NA.", call.=FALSE)
      } else {
        stop("Argument '", name, "' must be a probability between 0 and 1.", call.=FALSE)
      }
    }
  }
}

stopifNegative <- function(arg, allowNA=FALSE, allowZero=TRUE) {
  name <- deparse(substitute(arg))
  if(allowNA && all(is.na(arg))) {  # An all-NA vector is logical, but ok.
    # do nothing
  } else {
    if(!allowNA && any(is.na(arg)))
      stop("Argument '", name, "' must not contain NA or NaN.", call.=FALSE)
    if(!is.numeric(arg))
      stop("Argument '", name, "' must be numeric.", call.=FALSE)
    if(allowZero) {
      if(any(arg < 0, na.rm=TRUE)) {
        if(allowNA) {
          stop("Argument '", name, "' must be non-negative, or NA.", call.=FALSE)
        } else {
          stop("Argument '", name, "' must be non-negative.", call.=FALSE)
        }
      }
    } else {
      if(any(arg <= 0, na.rm=TRUE)) {
        if(allowNA) {
          stop("Argument '", name, "' must be greater than 0, or NA.", call.=FALSE)
        } else {
          stop("Argument '", name, "' must be greater than 0.", call.=FALSE)
        }
      }
    }
  }
}

# https://github.com/kenkellner/IPMbook/blob/main/R/stopifnot.R

# https://github.com/kenkellner/IPMbook/blob/main/R/BetaDist.R

getBeta2Par <- function(mean, sd) {
  stopifnotProbability(mean, allowNA=FALSE)
  stopifNegative(sd, allowNA=FALSE, allowZero=FALSE)

  nu <- mean * (1-mean) / sd^2 - 1
  if(any(nu <= 0)) {
    warning("sd is too large; some shape parameters will be NA.", call.=FALSE)
    nu[nu <= 0] <- NA
  }
  alpha <- mean * nu
  beta <- (1-mean) * nu
  cbind(shape1=alpha, shape2=beta)
}

dbeta2 <- function(x, mean, sd) {
  shapes <- getBeta2Par(mean,sd)
  return(dbeta(x, shapes[,1], shapes[,2]))
}

pbeta2 <- function(q, mean, sd, lower.tail=TRUE, log.p=FALSE) {
  shapes <- getBeta2Par(mean,sd)
  return(pbeta(q, shapes[,1], shapes[,2], lower.tail=lower.tail, log.p=log.p))
}
qbeta2 <- function(p, mean, sd, lower.tail=TRUE, log.p=FALSE) {
  shapes <- getBeta2Par(mean,sd)
  return(qbeta(p, shapes[,1], shapes[,2], lower.tail=lower.tail, log.p=log.p))
}

rbeta2 <- function(n, mean, sd) {
  shapes <- getBeta2Par(mean,sd)
  return(rbeta(n, shapes[,1], shapes[,2]))
}

# https://github.com/kenkellner/IPMbook/blob/main/R/BetaDist.R

Actual simulations

# https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R

# 3.3 Classical analysis of a matrix population model
# ===================================================

# 3.3.5 Analysis of a matrix population model with different sources of
#       stochasticity and parameter uncertainty
# ---------------------------------------------------------------------

# Define mean, measurement error and temporal variability of the demographic parameters
mean.sj <- 0.1163615          # Mean value of juv. survival
sd.sj.e <- 0.0554588        # Uncertainty of mean juv. survival as SD on natural scale
sd.sj.t <- plogis(0.03304038)         # population variability of juv. survival as SD on logit scale # Andrew, this might have to be plogis(0.03304038), but not sure, check the text.
mean.sa <- 0.5222         # Mean value of ad. survival
sd.sa.e <- 0.05769693        # Uncertainty of mean ad. survival as SD on natural scale
sd.sa.t <- plogis(0.0386221)         # Temporal variability of ad. survival as SD on logit scale
mean.f1 <- 1.4035          # Mean value of productivity of 1y females
sd.f1.e <- 0.1326032         # Uncertainty of mean productivity as SD on natural scale
sd.f1.t <- 0.1014479          # Temporal variability of productivity as SD on natural scale
mean.fa <- 1.4035          # Mean value of productivity of adult females
sd.fa.e <- 0.1326032         # Uncertainty of mean productivity as SD on natural scale
sd.fa.t <- 0.101          # Temporal variability of productivity as SD on natural scale
mean.imm <- 0.192587      # Mean value of adult immigration  
sd.imm.e <- 0.06674809    # Uncertainty of mean adult immigration as SD on natural scale
sd.imm.t <- 0.02872281    # population variability of adult immigration as SD on natrual scale

# Define the number of years with predictions and the Monte Carlo setting
T <- 10                # Number of years (projection time frame)
nsim <- 10000          # Number of replicate populations simulated

# Define population matrix and initial stage-specific population sizes
N <- array(NA, dim=c(2, T+1, nsim))
N[,1,] <- c(3,1)
r <- matrix(NA, nrow=T, ncol=nsim)
alive <- matrix(NA, nrow=T, ncol=nsim)
mean.r <- numeric(nsim)

# Project population
for (s in 1:nsim){ # Loop over replicate populations
  #if(s %% 100 == 0) {cat(paste("*** Simrep", s, "***\n")) }   # Counter
  # Generate a mean of the demographic rates (subject to measurement error)
  msj <- rbeta2(1, mean.sj, sd.sj.e)
  msa <- rbeta2(1, mean.sa, sd.sa.e)
  mf1 <- rnorm(1, mean.f1, sd.f1.e)
  mfa <- rnorm(1, mean.fa, sd.fa.e)
  mia <- rnorm(1, mean.imm, sd.imm.e)

  # Generate annual demographic rates (subject to temporal variability)
  sj <- plogis(rnorm(T, qlogis(msj), sd.sj.t))
  sa <- plogis(rnorm(T, qlogis(msa), sd.sa.t))
  f1 <- pmax(0, rnorm(T, mf1, sd.f1.t))                       # Avoids negative values
  fa <- pmax(0, rnorm(T, mfa, sd.fa.t))                       # Ditto
  ia <- pmax(0, rnorm(T, mia, sd.imm.t))                       # Ditto

  # Project population (include demographic stochasticity)
  for (t in 1:T){                                             # Loop over years
    N[1,t+1,s] <- rpois(1, sj[t] * (f1[t] * N[1,t,s] + fa[t] * N[2,t,s]))
    N[2,t+1,s] <- rbinom(1, sum(N[,t,s]), sa[t])+rpois(1,ia)
    if (sum(N[,t+1,s]) == 0) break
    #r[t,s] <- log(sum(N[,t+1,s])) - log(sum(N[,t,s]))
    alive[t,s] <- t
  } #t
  #mean.r[s] <- mean(r[min(alive[,s], na.rm=TRUE):max(alive[,s], na.rm=TRUE),s])
} 

#mean(mean.r)

#sd(mean.r)


not.extinct <- which(!is.na(alive[T,]))
mean(mean.r[not.extinct])

## [1] 0

sd(mean.r[not.extinct])

## [1] 0

# Extinction probability (after T years)
sum(is.na(alive[T,])) / nsim

## [1] 0.8031

rowSums(is.na(alive)) / nsim

##  [1] 0.0342 0.1379 0.2686 0.3941 0.5065 0.5914 0.6652 0.7232 0.7672 0.8031

df <- data.frame(
  totalcount=c(colSums(N[,,1:nsim])),
  time=rep(c(0:T),nsim),
  simnumber=gl(nsim, T+1)
)


# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# boxplot(c(colSums(N[,,1:nsim]))~factor(rep(c(0:T),nsim)))
 bxpobj<-boxplot(c(colSums(N[,,1:nsim]))~factor(rep(c(0:T),nsim)), plot=F)
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# points(bxpobj$stats[3,]~c(0:T), type="l", col="red", lwd=2)
# 
# par(mai=c(1,1,0.1,1))
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# points(bxpobj$stats[3,]~c(0:T), type="b", col="blue", lwd=2, pch=16)
# par(new=T)
# plot(c(rowSums(is.na(alive))/nsim)~c(1:T), type="b", col="red", bty="n", xlab="", ylab="", axes=F, ylim=c(0,1), xlim=c(0,T), lwd=2, pch=16)
# axis(side=4, at=pretty(range(c(0,1))))
# mtext("probability of extinction", side=4, line=3)

par(mai=c(1,1,0.1,1))
plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3, xlab="years of simulation", ylab="estimated population size")
points(bxpobj$stats[3,]~c(0:T), type="b", col="blue", lwd=2, pch=16)
par(new=T)
plot(c(rowSums(is.na(alive))/nsim)~c(1:T), type="b", col="red", bty="n", xlab="", ylab="", axes=F, ylim=c(0,1), xlim=c(0,T), lwd=2, pch=16)
axis(side=4, at=pretty(range(c(0,1))))
mtext("proportion of simulations reached N=0", side=4, line=3, col="red")

ggplot()+
 geom_line(data=df, aes(x=time, y=totalcount, group=simnumber), linewidth=1, color="gray", alpha=0.8) +
 geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
 #geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
 geom_label(aes(x=c(1:T), y=0, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
 scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
 theme_classic(base_size=22) +
 labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")

## Warning: Removed 40883 rows containing missing values or values outside the scale range
## (`geom_line()`).

ggplot()+
  geom_jitter(data=df, aes(x=time, y=totalcount, group=simnumber), shape=".", width=0.3, height=0.3) +
  geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
  #geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
  geom_label(aes(x=c(1:T), y=-0.5, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
  scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
  theme_classic(base_size=22)+
  labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")

## Warning: Removed 40883 rows containing missing values or values outside the scale range
## (`geom_point()`).

# ggplot()+
#   geom_hex(data=df, aes(x=time, y=totalcount, group=simnumber), bins=10) +
#   geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
#   #geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
#   geom_label(aes(x=c(1:T), y=-0.5, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
#   scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
#   theme_classic(base_size=22)+
#   labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")


# ~~~~ Plot graph with population growth rates (Fig. 3.15)  ~~~~
# op <- par(mar=c(4, 4, 2, 1), las=1, cex=1.1, "mfrow")
# layout(matrix(c(1, 1, 2, 3), 2, 2, byrow=TRUE), widths=c(1, 1), heights=c(1, 1), TRUE)
# plot(r[,1], type="l", lwd=0.5, ylab="Annual population growth rate", xlab="Time",
#     ylim=range(r[which(!is.na(alive))]), col="lightgrey", axes=FALSE)
# axis(1); axis(2)
# for (s in 2:nsim){
#   lines(r[!is.na(alive[,s]),s], lwd=0.5, col="lightgrey")
# }
# lines(apply(r, 1, mean, na.rm=TRUE), lwd=1.5)
# mtext("A", at=0, cex=1.5)
# a <- hist(mean.r, nclass=50, col="dodgerblue", main="",
#     xlab="Population growth rate", axes=FALSE) #xlim=c(-2.5, 0.1), 
# axis(1)
# axis(2, at=c(0, 5000, 10000, 15000, 20000, 25000, 30000), labels=c(0, 5, 10, 15, 20, 25, 30))
# mtext("B", at=a$mids[1], cex=1.5)
# a <- hist(mean.r[not.extinct], nclass=25, col="dodgerblue", main="",
#     xlab="Population growth rate", axes=FALSE)
# axis(1)
# axis(2, at=c(0, 2000, 4000, 6000, 8000, 10000), labels=c(0, 2, 4, 6, 8, 10))
# mtext("C", at=a$mids[1], cex=1.5)
# par(op)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R

Control Treatment simulations

# https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R

# 3.3 Classical analysis of a matrix population model
# ===================================================

# 3.3.5 Analysis of a matrix population model with different sources of
#       stochasticity and parameter uncertainty
# ---------------------------------------------------------------------

# Define mean, measurement error and temporal variability of the demographic parameters
mean.sj <- 0.1163615          # Mean value of juv. survival
sd.sj.e <- 0.0554588        # Uncertainty of mean juv. survival as SD on natural scale
sd.sj.t <- plogis(0.03304038)         # population variability of juv. survival as SD on logit scale # Andrew, this might have to be plogis(0.03304038), but not sure, check the text.
mean.sa <- 0.5222         # Mean value of ad. survival
sd.sa.e <- 0.05769693        # Uncertainty of mean ad. survival as SD on natural scale
sd.sa.t <- plogis(0.0386221)         # Temporal variability of ad. survival as SD on logit scale
mean.f1 <- 2.575          # Mean value of productivity of 1y females
sd.f1.e <- 0.1326032         # Uncertainty of mean productivity as SD on natural scale
sd.f1.t <- 0.1014479          # Temporal variability of productivity as SD on natural scale
mean.fa <- 2.575          # Mean value of productivity of adult females
sd.fa.e <- 0.1326032         # Uncertainty of mean productivity as SD on natural scale
sd.fa.t <- 0.101          # Temporal variability of productivity as SD on natural scale
mean.imm <- 0.192587      # Mean value of adult immigration  
sd.imm.e <- 0.06674809    # Uncertainty of mean adult immigration as SD on natural scale
sd.imm.t <- 0.02872281    # population variability of adult immigration as SD on natrual scale

# Define the number of years with predictions and the Monte Carlo setting
T <- 10                # Number of years (projection time frame)
nsim <- 10000          # Number of replicate populations simulated

# Define population matrix and initial stage-specific population sizes
N <- array(NA, dim=c(2, T+1, nsim))
N[,1,] <- c(3,1)
r <- matrix(NA, nrow=T, ncol=nsim)
alive <- matrix(NA, nrow=T, ncol=nsim)
mean.r <- numeric(nsim)

# Project population
for (s in 1:nsim){ # Loop over replicate populations
  #if(s %% 100 == 0) {cat(paste("*** Simrep", s, "***\n")) }   # Counter
  # Generate a mean of the demographic rates (subject to measurement error)
  msj <- rbeta2(1, mean.sj, sd.sj.e)
  msa <- rbeta2(1, mean.sa, sd.sa.e)
  mf1 <- rnorm(1, mean.f1, sd.f1.e)
  mfa <- rnorm(1, mean.fa, sd.fa.e)
  mia <- rnorm(1, mean.imm, sd.imm.e)

  # Generate annual demographic rates (subject to temporal variability)
  sj <- plogis(rnorm(T, qlogis(msj), sd.sj.t))
  sa <- plogis(rnorm(T, qlogis(msa), sd.sa.t))
  f1 <- pmax(0, rnorm(T, mf1, sd.f1.t))                       # Avoids negative values
  fa <- pmax(0, rnorm(T, mfa, sd.fa.t))                       # Ditto
  ia <- pmax(0, rnorm(T, mia, sd.imm.t))                       # Ditto

  # Project population (include demographic stochasticity)
  for (t in 1:T){                                             # Loop over years
    N[1,t+1,s] <- rpois(1, sj[t] * (f1[t] * N[1,t,s] + fa[t] * N[2,t,s]))
    N[2,t+1,s] <- rbinom(1, sum(N[,t,s]), sa[t])+rpois(1,ia)
    if (sum(N[,t+1,s]) == 0) break
    #r[t,s] <- log(sum(N[,t+1,s])) - log(sum(N[,t,s]))
    alive[t,s] <- t
  } #t
  #mean.r[s] <- mean(r[min(alive[,s], na.rm=TRUE):max(alive[,s], na.rm=TRUE),s])
} 

#mean(mean.r)

#sd(mean.r)


not.extinct <- which(!is.na(alive[T,]))
mean(mean.r[not.extinct])

## [1] 0

sd(mean.r[not.extinct])

## [1] 0

# Extinction probability (after T years)
sum(is.na(alive[T,])) / nsim

## [1] 0.5943

rowSums(is.na(alive)) / nsim

##  [1] 0.0254 0.0902 0.1752 0.2637 0.3443 0.4097 0.4697 0.5164 0.5589 0.5943

df <- data.frame(
  totalcount=c(colSums(N[,,1:nsim])),
  time=rep(c(0:T),nsim),
  simnumber=gl(nsim, T+1)
)


# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# boxplot(c(colSums(N[,,1:nsim]))~factor(rep(c(0:T),nsim)))
 bxpobj<-boxplot(c(colSums(N[,,1:nsim]))~factor(rep(c(0:T),nsim)), plot=F)
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# points(bxpobj$stats[3,]~c(0:T), type="l", col="red", lwd=2)
# 
# par(mai=c(1,1,0.1,1))
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# points(bxpobj$stats[3,]~c(0:T), type="b", col="blue", lwd=2, pch=16)
# par(new=T)
# plot(c(rowSums(is.na(alive))/nsim)~c(1:T), type="b", col="red", bty="n", xlab="", ylab="", axes=F, ylim=c(0,1), xlim=c(0,T), lwd=2, pch=16)
# axis(side=4, at=pretty(range(c(0,1))))
# mtext("probability of extinction", side=4, line=3)

par(mai=c(1,1,0.1,1))
plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3, xlab="years of simulation", ylab="estimated population size")
points(bxpobj$stats[3,]~c(0:T), type="b", col="blue", lwd=2, pch=16)
par(new=T)
plot(c(rowSums(is.na(alive))/nsim)~c(1:T), type="b", col="red", bty="n", xlab="", ylab="", axes=F, ylim=c(0,1), xlim=c(0,T), lwd=2, pch=16)
axis(side=4, at=pretty(range(c(0,1))))
mtext("proportion of simulations reached N=0", side=4, line=3, col="red")

ggplot()+
 geom_line(data=df, aes(x=time, y=totalcount, group=simnumber), linewidth=1, color="gray", alpha=0.8) +
 geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
 #geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
 geom_label(aes(x=c(1:T), y=0, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
 scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
 theme_classic(base_size=22) +
 labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")

## Warning: Removed 28535 rows containing missing values or values outside the scale range
## (`geom_line()`).

ggplot()+
  geom_jitter(data=df, aes(x=time, y=totalcount, group=simnumber), shape=".", width=0.3, height=0.3) +
  geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
  #geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
  geom_label(aes(x=c(1:T), y=-0.5, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
  scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
  theme_classic(base_size=22)+
  labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")

## Warning: Removed 28535 rows containing missing values or values outside the scale range
## (`geom_point()`).

# ggplot()+
#   geom_hex(data=df, aes(x=time, y=totalcount, group=simnumber), bins=10) +
#   geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
#   #geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
#   geom_label(aes(x=c(1:T), y=-0.5, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
#   scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
#   theme_classic(base_size=22)+
#   labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")


# ~~~~ Plot graph with population growth rates (Fig. 3.15)  ~~~~
# op <- par(mar=c(4, 4, 2, 1), las=1, cex=1.1, "mfrow")
# layout(matrix(c(1, 1, 2, 3), 2, 2, byrow=TRUE), widths=c(1, 1), heights=c(1, 1), TRUE)
# plot(r[,1], type="l", lwd=0.5, ylab="Annual population growth rate", xlab="Time",
#     ylim=range(r[which(!is.na(alive))]), col="lightgrey", axes=FALSE)
# axis(1); axis(2)
# for (s in 2:nsim){
#   lines(r[!is.na(alive[,s]),s], lwd=0.5, col="lightgrey")
# }
# lines(apply(r, 1, mean, na.rm=TRUE), lwd=1.5)
# mtext("A", at=0, cex=1.5)
# a <- hist(mean.r, nclass=50, col="dodgerblue", main="",
#     xlab="Population growth rate", axes=FALSE) #xlim=c(-2.5, 0.1), 
# axis(1)
# axis(2, at=c(0, 5000, 10000, 15000, 20000, 25000, 30000), labels=c(0, 5, 10, 15, 20, 25, 30))
# mtext("B", at=a$mids[1], cex=1.5)
# a <- hist(mean.r[not.extinct], nclass=25, col="dodgerblue", main="",
#     xlab="Population growth rate", axes=FALSE)
# axis(1)
# axis(2, at=c(0, 2000, 4000, 6000, 8000, 10000), labels=c(0, 2, 4, 6, 8, 10))
# mtext("C", at=a$mids[1], cex=1.5)
# par(op)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R

Low Treatment simulations

# https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R

# 3.3 Classical analysis of a matrix population model
# ===================================================

# 3.3.5 Analysis of a matrix population model with different sources of
#       stochasticity and parameter uncertainty
# ---------------------------------------------------------------------

# Define mean, measurement error and temporal variability of the demographic parameters
mean.sj <- 0.1163615          # Mean value of juv. survival
sd.sj.e <- 0.0554588        # Uncertainty of mean juv. survival as SD on natural scale
sd.sj.t <- plogis(0.03304038)         # population variability of juv. survival as SD on logit scale # Andrew, this might have to be plogis(0.03304038), but not sure, check the text.
mean.sa <- 0.5222         # Mean value of ad. survival
sd.sa.e <- 0.05769693        # Uncertainty of mean ad. survival as SD on natural scale
sd.sa.t <- plogis(0.0386221)         # Temporal variability of ad. survival as SD on logit scale
mean.f1 <- 2.264          # Mean value of productivity of 1y females
sd.f1.e <- 0.1326032         # Uncertainty of mean productivity as SD on natural scale
sd.f1.t <- 0.1014479          # Temporal variability of productivity as SD on natural scale
mean.fa <- 2.264          # Mean value of productivity of adult females
sd.fa.e <- 0.1326032         # Uncertainty of mean productivity as SD on natural scale
sd.fa.t <- 0.101          # Temporal variability of productivity as SD on natural scale
mean.imm <- 0.192587      # Mean value of adult immigration  
sd.imm.e <- 0.06674809    # Uncertainty of mean adult immigration as SD on natural scale
sd.imm.t <- 0.02872281    # population variability of adult immigration as SD on natrual scale

# Define the number of years with predictions and the Monte Carlo setting
T <- 10                # Number of years (projection time frame)
nsim <- 10000          # Number of replicate populations simulated

# Define population matrix and initial stage-specific population sizes
N <- array(NA, dim=c(2, T+1, nsim))
N[,1,] <- c(3,1)
r <- matrix(NA, nrow=T, ncol=nsim)
alive <- matrix(NA, nrow=T, ncol=nsim)
mean.r <- numeric(nsim)

# Project population
for (s in 1:nsim){ # Loop over replicate populations
  #if(s %% 100 == 0) {cat(paste("*** Simrep", s, "***\n")) }   # Counter
  # Generate a mean of the demographic rates (subject to measurement error)
  msj <- rbeta2(1, mean.sj, sd.sj.e)
  msa <- rbeta2(1, mean.sa, sd.sa.e)
  mf1 <- rnorm(1, mean.f1, sd.f1.e)
  mfa <- rnorm(1, mean.fa, sd.fa.e)
  mia <- rnorm(1, mean.imm, sd.imm.e)

  # Generate annual demographic rates (subject to temporal variability)
  sj <- plogis(rnorm(T, qlogis(msj), sd.sj.t))
  sa <- plogis(rnorm(T, qlogis(msa), sd.sa.t))
  f1 <- pmax(0, rnorm(T, mf1, sd.f1.t))                       # Avoids negative values
  fa <- pmax(0, rnorm(T, mfa, sd.fa.t))                       # Ditto
  ia <- pmax(0, rnorm(T, mia, sd.imm.t))                       # Ditto

  # Project population (include demographic stochasticity)
  for (t in 1:T){                                             # Loop over years
    N[1,t+1,s] <- rpois(1, sj[t] * (f1[t] * N[1,t,s] + fa[t] * N[2,t,s]))
    N[2,t+1,s] <- rbinom(1, sum(N[,t,s]), sa[t])+rpois(1,ia)
    if (sum(N[,t+1,s]) == 0) break
    #r[t,s] <- log(sum(N[,t+1,s])) - log(sum(N[,t,s]))
    alive[t,s] <- t
  } #t
  #mean.r[s] <- mean(r[min(alive[,s], na.rm=TRUE):max(alive[,s], na.rm=TRUE),s])
} 

#mean(mean.r)

#sd(mean.r)


not.extinct <- which(!is.na(alive[T,]))
mean(mean.r[not.extinct])

## [1] 0

sd(mean.r[not.extinct])

## [1] 0

# Extinction probability (after T years)
sum(is.na(alive[T,])) / nsim

## [1] 0.6536

rowSums(is.na(alive)) / nsim

##  [1] 0.0270 0.1089 0.2069 0.2994 0.3825 0.4560 0.5194 0.5701 0.6143 0.6536

df <- data.frame(
  totalcount=c(colSums(N[,,1:nsim])),
  time=rep(c(0:T),nsim),
  simnumber=gl(nsim, T+1)
)


# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# boxplot(c(colSums(N[,,1:nsim]))~factor(rep(c(0:T),nsim)))
 bxpobj<-boxplot(c(colSums(N[,,1:nsim]))~factor(rep(c(0:T),nsim)), plot=F)
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# points(bxpobj$stats[3,]~c(0:T), type="l", col="red", lwd=2)
# 
# par(mai=c(1,1,0.1,1))
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# points(bxpobj$stats[3,]~c(0:T), type="b", col="blue", lwd=2, pch=16)
# par(new=T)
# plot(c(rowSums(is.na(alive))/nsim)~c(1:T), type="b", col="red", bty="n", xlab="", ylab="", axes=F, ylim=c(0,1), xlim=c(0,T), lwd=2, pch=16)
# axis(side=4, at=pretty(range(c(0,1))))
# mtext("probability of extinction", side=4, line=3)

par(mai=c(1,1,0.1,1))
plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3, xlab="years of simulation", ylab="estimated population size")
points(bxpobj$stats[3,]~c(0:T), type="b", col="blue", lwd=2, pch=16)
par(new=T)
plot(c(rowSums(is.na(alive))/nsim)~c(1:T), type="b", col="red", bty="n", xlab="", ylab="", axes=F, ylim=c(0,1), xlim=c(0,T), lwd=2, pch=16)
axis(side=4, at=pretty(range(c(0,1))))
mtext("proportion of simulations reached N=0", side=4, line=3, col="red")

ggplot()+
 geom_line(data=df, aes(x=time, y=totalcount, group=simnumber), linewidth=1, color="gray", alpha=0.8) +
 geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
 #geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
 geom_label(aes(x=c(1:T), y=0, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
 scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
 theme_classic(base_size=22) +
 labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")

## Warning: Removed 31845 rows containing missing values or values outside the scale range
## (`geom_line()`).

ggplot()+
  geom_jitter(data=df, aes(x=time, y=totalcount, group=simnumber), shape=".", width=0.3, height=0.3) +
  geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
  #geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
  geom_label(aes(x=c(1:T), y=-0.5, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
  scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
  theme_classic(base_size=22)+
  labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")

## Warning: Removed 31845 rows containing missing values or values outside the scale range
## (`geom_point()`).

# ggplot()+
#   geom_hex(data=df, aes(x=time, y=totalcount, group=simnumber), bins=10) +
#   geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
#   #geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
#   geom_label(aes(x=c(1:T), y=-0.5, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
#   scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
#   theme_classic(base_size=22)+
#   labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")


# ~~~~ Plot graph with population growth rates (Fig. 3.15)  ~~~~
# op <- par(mar=c(4, 4, 2, 1), las=1, cex=1.1, "mfrow")
# layout(matrix(c(1, 1, 2, 3), 2, 2, byrow=TRUE), widths=c(1, 1), heights=c(1, 1), TRUE)
# plot(r[,1], type="l", lwd=0.5, ylab="Annual population growth rate", xlab="Time",
#     ylim=range(r[which(!is.na(alive))]), col="lightgrey", axes=FALSE)
# axis(1); axis(2)
# for (s in 2:nsim){
#   lines(r[!is.na(alive[,s]),s], lwd=0.5, col="lightgrey")
# }
# lines(apply(r, 1, mean, na.rm=TRUE), lwd=1.5)
# mtext("A", at=0, cex=1.5)
# a <- hist(mean.r, nclass=50, col="dodgerblue", main="",
#     xlab="Population growth rate", axes=FALSE) #xlim=c(-2.5, 0.1), 
# axis(1)
# axis(2, at=c(0, 5000, 10000, 15000, 20000, 25000, 30000), labels=c(0, 5, 10, 15, 20, 25, 30))
# mtext("B", at=a$mids[1], cex=1.5)
# a <- hist(mean.r[not.extinct], nclass=25, col="dodgerblue", main="",
#     xlab="Population growth rate", axes=FALSE)
# axis(1)
# axis(2, at=c(0, 2000, 4000, 6000, 8000, 10000), labels=c(0, 2, 4, 6, 8, 10))
# mtext("C", at=a$mids[1], cex=1.5)
# par(op)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R

Medium Treatment simulations

# https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R

# 3.3 Classical analysis of a matrix population model
# ===================================================

# 3.3.5 Analysis of a matrix population model with different sources of
#       stochasticity and parameter uncertainty
# ---------------------------------------------------------------------

# Define mean, measurement error and temporal variability of the demographic parameters
mean.sj <- 0.1163615          # Mean value of juv. survival
sd.sj.e <- 0.0554588        # Uncertainty of mean juv. survival as SD on natural scale
sd.sj.t <- plogis(0.03304038)         # population variability of juv. survival as SD on logit scale # Andrew, this might have to be plogis(0.03304038), but not sure, check the text.
mean.sa <- 0.5222         # Mean value of ad. survival
sd.sa.e <- 0.05769693        # Uncertainty of mean ad. survival as SD on natural scale
sd.sa.t <- plogis(0.0386221)         # Temporal variability of ad. survival as SD on logit scale
mean.f1 <- 1.605          # Mean value of productivity of 1y females
sd.f1.e <- 0.1326032         # Uncertainty of mean productivity as SD on natural scale
sd.f1.t <- 0.1014479          # Temporal variability of productivity as SD on natural scale
mean.fa <- 1.605          # Mean value of productivity of adult females
sd.fa.e <- 0.1326032         # Uncertainty of mean productivity as SD on natural scale
sd.fa.t <- 0.101          # Temporal variability of productivity as SD on natural scale
mean.imm <- 0.192587      # Mean value of adult immigration  
sd.imm.e <- 0.06674809    # Uncertainty of mean adult immigration as SD on natural scale
sd.imm.t <- 0.02872281    # population variability of adult immigration as SD on natrual scale

# Define the number of years with predictions and the Monte Carlo setting
T <- 10                # Number of years (projection time frame)
nsim <- 10000          # Number of replicate populations simulated

# Define population matrix and initial stage-specific population sizes
N <- array(NA, dim=c(2, T+1, nsim))
N[,1,] <- c(3,1)
r <- matrix(NA, nrow=T, ncol=nsim)
alive <- matrix(NA, nrow=T, ncol=nsim)
mean.r <- numeric(nsim)

# Project population
for (s in 1:nsim){ # Loop over replicate populations
  #if(s %% 100 == 0) {cat(paste("*** Simrep", s, "***\n")) }   # Counter
  # Generate a mean of the demographic rates (subject to measurement error)
  msj <- rbeta2(1, mean.sj, sd.sj.e)
  msa <- rbeta2(1, mean.sa, sd.sa.e)
  mf1 <- rnorm(1, mean.f1, sd.f1.e)
  mfa <- rnorm(1, mean.fa, sd.fa.e)
  mia <- rnorm(1, mean.imm, sd.imm.e)

  # Generate annual demographic rates (subject to temporal variability)
  sj <- plogis(rnorm(T, qlogis(msj), sd.sj.t))
  sa <- plogis(rnorm(T, qlogis(msa), sd.sa.t))
  f1 <- pmax(0, rnorm(T, mf1, sd.f1.t))                       # Avoids negative values
  fa <- pmax(0, rnorm(T, mfa, sd.fa.t))                       # Ditto
  ia <- pmax(0, rnorm(T, mia, sd.imm.t))                       # Ditto

  # Project population (include demographic stochasticity)
  for (t in 1:T){                                             # Loop over years
    N[1,t+1,s] <- rpois(1, sj[t] * (f1[t] * N[1,t,s] + fa[t] * N[2,t,s]))
    N[2,t+1,s] <- rbinom(1, sum(N[,t,s]), sa[t])+rpois(1,ia)
    if (sum(N[,t+1,s]) == 0) break
    #r[t,s] <- log(sum(N[,t+1,s])) - log(sum(N[,t,s]))
    alive[t,s] <- t
  } #t
  #mean.r[s] <- mean(r[min(alive[,s], na.rm=TRUE):max(alive[,s], na.rm=TRUE),s])
} 

#mean(mean.r)

#sd(mean.r)


not.extinct <- which(!is.na(alive[T,]))
mean(mean.r[not.extinct])

## [1] 0

sd(mean.r[not.extinct])

## [1] 0

# Extinction probability (after T years)
sum(is.na(alive[T,])) / nsim

## [1] 0.7736

rowSums(is.na(alive)) / nsim

##  [1] 0.0300 0.1266 0.2507 0.3666 0.4616 0.5496 0.6284 0.6869 0.7358 0.7736

df <- data.frame(
  totalcount=c(colSums(N[,,1:nsim])),
  time=rep(c(0:T),nsim),
  simnumber=gl(nsim, T+1)
)


# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# boxplot(c(colSums(N[,,1:nsim]))~factor(rep(c(0:T),nsim)))
 bxpobj<-boxplot(c(colSums(N[,,1:nsim]))~factor(rep(c(0:T),nsim)), plot=F)
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# points(bxpobj$stats[3,]~c(0:T), type="l", col="red", lwd=2)
# 
# par(mai=c(1,1,0.1,1))
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# points(bxpobj$stats[3,]~c(0:T), type="b", col="blue", lwd=2, pch=16)
# par(new=T)
# plot(c(rowSums(is.na(alive))/nsim)~c(1:T), type="b", col="red", bty="n", xlab="", ylab="", axes=F, ylim=c(0,1), xlim=c(0,T), lwd=2, pch=16)
# axis(side=4, at=pretty(range(c(0,1))))
# mtext("probability of extinction", side=4, line=3)

par(mai=c(1,1,0.1,1))
plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3, xlab="years of simulation", ylab="estimated population size")
points(bxpobj$stats[3,]~c(0:T), type="b", col="blue", lwd=2, pch=16)
par(new=T)
plot(c(rowSums(is.na(alive))/nsim)~c(1:T), type="b", col="red", bty="n", xlab="", ylab="", axes=F, ylim=c(0,1), xlim=c(0,T), lwd=2, pch=16)
axis(side=4, at=pretty(range(c(0,1))))
mtext("proportion of simulations reached N=0", side=4, line=3, col="red")

ggplot()+
 geom_line(data=df, aes(x=time, y=totalcount, group=simnumber), linewidth=1, color="gray", alpha=0.8) +
 geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
 #geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
 geom_label(aes(x=c(1:T), y=0, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
 scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
 theme_classic(base_size=22) +
 labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")

## Warning: Removed 38362 rows containing missing values or values outside the scale range
## (`geom_line()`).

ggplot()+
  geom_jitter(data=df, aes(x=time, y=totalcount, group=simnumber), shape=".", width=0.3, height=0.3) +
  geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
  #geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
  geom_label(aes(x=c(1:T), y=-0.5, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
  scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
  theme_classic(base_size=22)+
  labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")

## Warning: Removed 38362 rows containing missing values or values outside the scale range
## (`geom_point()`).

# ggplot()+
#   geom_hex(data=df, aes(x=time, y=totalcount, group=simnumber), bins=10) +
#   geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
#   #geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
#   geom_label(aes(x=c(1:T), y=-0.5, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
#   scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
#   theme_classic(base_size=22)+
#   labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")


# ~~~~ Plot graph with population growth rates (Fig. 3.15)  ~~~~
# op <- par(mar=c(4, 4, 2, 1), las=1, cex=1.1, "mfrow")
# layout(matrix(c(1, 1, 2, 3), 2, 2, byrow=TRUE), widths=c(1, 1), heights=c(1, 1), TRUE)
# plot(r[,1], type="l", lwd=0.5, ylab="Annual population growth rate", xlab="Time",
#     ylim=range(r[which(!is.na(alive))]), col="lightgrey", axes=FALSE)
# axis(1); axis(2)
# for (s in 2:nsim){
#   lines(r[!is.na(alive[,s]),s], lwd=0.5, col="lightgrey")
# }
# lines(apply(r, 1, mean, na.rm=TRUE), lwd=1.5)
# mtext("A", at=0, cex=1.5)
# a <- hist(mean.r, nclass=50, col="dodgerblue", main="",
#     xlab="Population growth rate", axes=FALSE) #xlim=c(-2.5, 0.1), 
# axis(1)
# axis(2, at=c(0, 5000, 10000, 15000, 20000, 25000, 30000), labels=c(0, 5, 10, 15, 20, 25, 30))
# mtext("B", at=a$mids[1], cex=1.5)
# a <- hist(mean.r[not.extinct], nclass=25, col="dodgerblue", main="",
#     xlab="Population growth rate", axes=FALSE)
# axis(1)
# axis(2, at=c(0, 2000, 4000, 6000, 8000, 10000), labels=c(0, 2, 4, 6, 8, 10))
# mtext("C", at=a$mids[1], cex=1.5)
# par(op)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R

High Treatment simulations

# https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R

# 3.3 Classical analysis of a matrix population model
# ===================================================

# 3.3.5 Analysis of a matrix population model with different sources of
#       stochasticity and parameter uncertainty
# ---------------------------------------------------------------------

# Define mean, measurement error and temporal variability of the demographic parameters
mean.sj <- 0.1163615          # Mean value of juv. survival
sd.sj.e <- 0.0554588        # Uncertainty of mean juv. survival as SD on natural scale
sd.sj.t <- plogis(0.03304038)         # population variability of juv. survival as SD on logit scale # Andrew, this might have to be plogis(0.03304038), but not sure, check the text.
mean.sa <- 0.5222         # Mean value of ad. survival
sd.sa.e <- 0.05769693        # Uncertainty of mean ad. survival as SD on natural scale
sd.sa.t <- plogis(0.0386221)         # Temporal variability of ad. survival as SD on logit scale
mean.f1 <- 0.343          # Mean value of productivity of 1y females
sd.f1.e <- 0.1326032         # Uncertainty of mean productivity as SD on natural scale
sd.f1.t <- 0.1014479          # Temporal variability of productivity as SD on natural scale
mean.fa <- 0.343          # Mean value of productivity of adult females
sd.fa.e <- 0.1326032         # Uncertainty of mean productivity as SD on natural scale
sd.fa.t <- 0.101          # Temporal variability of productivity as SD on natural scale
mean.imm <- 0.192587      # Mean value of adult immigration  
sd.imm.e <- 0.06674809    # Uncertainty of mean adult immigration as SD on natural scale
sd.imm.t <- 0.02872281    # population variability of adult immigration as SD on natrual scale

# Define the number of years with predictions and the Monte Carlo setting
T <- 10                # Number of years (projection time frame)
nsim <- 10000          # Number of replicate populations simulated

# Define population matrix and initial stage-specific population sizes
N <- array(NA, dim=c(2, T+1, nsim))
N[,1,] <- c(3,1)
r <- matrix(NA, nrow=T, ncol=nsim)
alive <- matrix(NA, nrow=T, ncol=nsim)
mean.r <- numeric(nsim)

# Project population
for (s in 1:nsim){ # Loop over replicate populations
  #if(s %% 100 == 0) {cat(paste("*** Simrep", s, "***\n")) }   # Counter
  # Generate a mean of the demographic rates (subject to measurement error)
  msj <- rbeta2(1, mean.sj, sd.sj.e)
  msa <- rbeta2(1, mean.sa, sd.sa.e)
  mf1 <- rnorm(1, mean.f1, sd.f1.e)
  mfa <- rnorm(1, mean.fa, sd.fa.e)
  mia <- rnorm(1, mean.imm, sd.imm.e)

  # Generate annual demographic rates (subject to temporal variability)
  sj <- plogis(rnorm(T, qlogis(msj), sd.sj.t))
  sa <- plogis(rnorm(T, qlogis(msa), sd.sa.t))
  f1 <- pmax(0, rnorm(T, mf1, sd.f1.t))                       # Avoids negative values
  fa <- pmax(0, rnorm(T, mfa, sd.fa.t))                       # Ditto
  ia <- pmax(0, rnorm(T, mia, sd.imm.t))                       # Ditto

  # Project population (include demographic stochasticity)
  for (t in 1:T){                                             # Loop over years
    N[1,t+1,s] <- rpois(1, sj[t] * (f1[t] * N[1,t,s] + fa[t] * N[2,t,s]))
    N[2,t+1,s] <- rbinom(1, sum(N[,t,s]), sa[t])+rpois(1,ia)
    if (sum(N[,t+1,s]) == 0) break
    #r[t,s] <- log(sum(N[,t+1,s])) - log(sum(N[,t,s]))
    alive[t,s] <- t
  } #t
  #mean.r[s] <- mean(r[min(alive[,s], na.rm=TRUE):max(alive[,s], na.rm=TRUE),s])
} 

#mean(mean.r)

#sd(mean.r)


not.extinct <- which(!is.na(alive[T,]))
mean(mean.r[not.extinct])

## [1] 0

sd(mean.r[not.extinct])

## [1] 0

# Extinction probability (after T years)
sum(is.na(alive[T,])) / nsim

## [1] 0.9372

rowSums(is.na(alive)) / nsim

##  [1] 0.0554 0.2181 0.4013 0.5635 0.6808 0.7702 0.8315 0.8801 0.9137 0.9372

df <- data.frame(
  totalcount=c(colSums(N[,,1:nsim])),
  time=rep(c(0:T),nsim),
  simnumber=gl(nsim, T+1)
)


# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# boxplot(c(colSums(N[,,1:nsim]))~factor(rep(c(0:T),nsim)))
 bxpobj<-boxplot(c(colSums(N[,,1:nsim]))~factor(rep(c(0:T),nsim)), plot=F)
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# points(bxpobj$stats[3,]~c(0:T), type="l", col="red", lwd=2)
# 
# par(mai=c(1,1,0.1,1))
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# points(bxpobj$stats[3,]~c(0:T), type="b", col="blue", lwd=2, pch=16)
# par(new=T)
# plot(c(rowSums(is.na(alive))/nsim)~c(1:T), type="b", col="red", bty="n", xlab="", ylab="", axes=F, ylim=c(0,1), xlim=c(0,T), lwd=2, pch=16)
# axis(side=4, at=pretty(range(c(0,1))))
# mtext("probability of extinction", side=4, line=3)

par(mai=c(1,1,0.1,1))
plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3, xlab="years of simulation", ylab="estimated population size")
points(bxpobj$stats[3,]~c(0:T), type="b", col="blue", lwd=2, pch=16)
par(new=T)
plot(c(rowSums(is.na(alive))/nsim)~c(1:T), type="b", col="red", bty="n", xlab="", ylab="", axes=F, ylim=c(0,1), xlim=c(0,T), lwd=2, pch=16)
axis(side=4, at=pretty(range(c(0,1))))
mtext("proportion of simulations reached N=0", side=4, line=3, col="red")

ggplot()+
 geom_line(data=df, aes(x=time, y=totalcount, group=simnumber), linewidth=1, color="gray", alpha=0.8) +
 geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
 #geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
 geom_label(aes(x=c(1:T), y=0, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
 scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
 theme_classic(base_size=22) +
 labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")

## Warning: Removed 53146 rows containing missing values or values outside the scale range
## (`geom_line()`).

ggplot()+
  geom_jitter(data=df, aes(x=time, y=totalcount, group=simnumber), shape=".", width=0.3, height=0.3) +
  geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
  #geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
  geom_label(aes(x=c(1:T), y=-0.5, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
  scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
  theme_classic(base_size=22)+
  labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")

## Warning: Removed 53146 rows containing missing values or values outside the scale range
## (`geom_point()`).

# ggplot()+
#   geom_hex(data=df, aes(x=time, y=totalcount, group=simnumber), bins=10) +
#   geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
#   #geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
#   geom_label(aes(x=c(1:T), y=-0.5, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
#   scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
#   theme_classic(base_size=22)+
#   labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")


# ~~~~ Plot graph with population growth rates (Fig. 3.15)  ~~~~
# op <- par(mar=c(4, 4, 2, 1), las=1, cex=1.1, "mfrow")
# layout(matrix(c(1, 1, 2, 3), 2, 2, byrow=TRUE), widths=c(1, 1), heights=c(1, 1), TRUE)
# plot(r[,1], type="l", lwd=0.5, ylab="Annual population growth rate", xlab="Time",
#     ylim=range(r[which(!is.na(alive))]), col="lightgrey", axes=FALSE)
# axis(1); axis(2)
# for (s in 2:nsim){
#   lines(r[!is.na(alive[,s]),s], lwd=0.5, col="lightgrey")
# }
# lines(apply(r, 1, mean, na.rm=TRUE), lwd=1.5)
# mtext("A", at=0, cex=1.5)
# a <- hist(mean.r, nclass=50, col="dodgerblue", main="",
#     xlab="Population growth rate", axes=FALSE) #xlim=c(-2.5, 0.1), 
# axis(1)
# axis(2, at=c(0, 5000, 10000, 15000, 20000, 25000, 30000), labels=c(0, 5, 10, 15, 20, 25, 30))
# mtext("B", at=a$mids[1], cex=1.5)
# a <- hist(mean.r[not.extinct], nclass=25, col="dodgerblue", main="",
#     xlab="Population growth rate", axes=FALSE)
# axis(1)
# axis(2, at=c(0, 2000, 4000, 6000, 8000, 10000), labels=c(0, 2, 4, 6, 8, 10))
# mtext("C", at=a$mids[1], cex=1.5)
# par(op)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R

treatment comparisons

# # for actual treatment fecundity values
# contm2 <- out1$mean$lambda[1] #2.575
# lowm2 <- out1$mean$lambda[7] #2.264
# medm2 <- out1$mean$lambda[14] #1.605
# highm2 <- out1$mean$lambda[20] #0.343
# alpha <- out1$mean$alpha #0.9422234
# beta <- out1$mean$beta #-1.073293

probofextinctdf <- data.frame(
  treatment=c(rep("Control",10),rep("Low",10),rep("Medium", 10), rep("High",10)),
  time=rep(1:10, 4),
  probextinct=c(0.0241,0.0896, 0.1795, 0.2620, 0.3443, 0.4128, 0.4762, 0.5264, 0.5677, 0.6027,
                0.0251, 0.1022, 0.2020, 0.2966, 0.3807, 0.4602, 0.5222, 0.5747, 0.6169, 0.6544,
                0.0331, 0.1322, 0.2553, 0.3724, 0.4682, 0.5516, 0.6191, 0.6773, 0.7268, 0.7654,
                0.0543, 0.2188, 0.4035, 0.5623, 0.6821, 0.7708, 0.8324, 0.8764, 0.9078, 0.9308),
  fecundity=c(rep(2.575, 10),rep(2.264,10),rep(1.605,10),rep(0.343,10))
)

ggplot()+
  geom_point(data=probofextinctdf, aes(x=time, y=probextinct, color=factor(treatment, levels=c("Control","Low","Medium","High"))))+
  geom_line(data=probofextinctdf, aes(x=time, y=probextinct, color=factor(treatment, levels=c("Control","Low","Medium","High"))))+
  labs(color="", x="Year of simulation", y="Proportion of simulations with N=0")+
  scale_y_continuous(breaks=c(0, 0.25, 0.5, 0.75, 1.0), limits=c(0,1))+
  scale_x_continuous(breaks=c(1:10))+
  theme_classic()

plot(probextinct~fecundity, data=probofextinctdf)
abline(lm(probextinct~fecundity, data=probofextinctdf))
abline(lm(probextinct~fecundity, data=probofextinctdf|>filter(time==10)), col="red")

summary(lm(probextinct~fecundity+time, data=probofextinctdf))

## 
## Call:
## lm(formula = probextinct ~ fecundity + time, data = probofextinctdf)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.20902 -0.02568  0.01506  0.03662  0.10939 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.226164   0.029695   7.616 4.35e-09 ***
## fecundity   -0.123788   0.011829 -10.465 1.31e-12 ***
## time         0.079617   0.003527  22.572  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06407 on 37 degrees of freedom
## Multiple R-squared:  0.9436, Adjusted R-squared:  0.9406 
## F-statistic: 309.5 on 2 and 37 DF,  p-value: < 2.2e-16

summary(lm(probextinct~fecundity, data=probofextinctdf|>filter(time==10)))

## 
## Call:
## lm(formula = probextinct ~ fecundity, data = filter(probofextinctdf, 
##     time == 10))
## 
## Residuals:
##          1          2          3          4 
## -0.0069874 -0.0008397  0.0136363 -0.0058092 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.986849   0.012874   76.66  0.00017 ***
## fecundity   -0.146470   0.006773  -21.62  0.00213 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0116 on 2 degrees of freedom
## Multiple R-squared:  0.9957, Adjusted R-squared:  0.9936 
## F-statistic: 467.6 on 1 and 2 DF,  p-value: 0.002132

summary(lm(probextinct~fecundity, data=probofextinctdf|>filter(time==5)))

## 
## Call:
## lm(formula = probextinct ~ fecundity, data = filter(probofextinctdf, 
##     time == 5))
## 
## Residuals:
##         1         2         3         4 
##  0.009575 -0.001511 -0.014634  0.006571 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.727902   0.014706   49.50 0.000408 ***
## fecundity   -0.152690   0.007737  -19.73 0.002558 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01325 on 2 degrees of freedom
## Multiple R-squared:  0.9949, Adjusted R-squared:  0.9923 
## F-statistic: 389.4 on 1 and 2 DF,  p-value: 0.002558

# -0.20*(1/-0.1526)
# = 1.310616
# A reduction in fecundity of 1.31 is expected to lead to a 0.2 increase in risk of N=0 in the 5th of 10 years

# 2.5750723-1.310616=1.264456

# alpha = 0.9422234
# beta = -1.073293

# (log(1.264456)-0.9422234)/-1.073293 =
# 0.6592621